%%%Main File for BCDR Simulations - gains from industrial policy only

%This file computes the gains from optimal industrial policy under the
%constraint of zero export taxes (Table 3, column 2 in BCDR)

%CES version

clear all;


K=32; %number of sectors
N=62; %includes ROW
K_man=15;    %number of manufacturing sectors

%%%Load data
%%%Data is organized as follows: first column origin country codes, second column
%%%sector codes, third column sectoral value added, fourth column aggregate origin trade balance,
%%%fifth column origin-sector expenditurs, 6th column starts the bilateral
%%%trade flows

%Note: data is sorted by sector, then origin.

A=importdata('simulation_ge_data.csv');
B=importdata('sim_params_no_int.csv'); %these are the estimated elasticities

%Data vectors
X_ink = A(:,6:5+N); %Matrix of trade flows, (NxK)xN
D_i=A(1:N,4); %vector of trade balances, Nx1
V_ik=A(:,3); %vector of sectoral value added, (NxK)x1
expenditures_ik=A(:,5); %vector of expenditures, (NxK)x1
countryid=A(1:N,1); %vector of country ids, Nx1

%Parameter vectors
theta_nm=4.32; %non-manufacturing sector TE set to median of manufacturing
h_theta = theta_nm*[ones(2,1);zeros(K_man,1);ones(K-2-K_man,1)]+[zeros(2,1);B(:,2);zeros(K-2-K_man,1)]; %trade elasticities
h_gamma = [zeros(2,1);B(:,1);zeros(K-2-K_man,1)]; %Heterogeneous scale elasticities, non-manufacturing set to zero
gammaXtheta=h_theta.*h_gamma;
rho=.87; %Upper tier elasticity of substitution

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Reshaping data vectors

%Reshape trade matrix to be NxNxK (3 dimensional array)
%rows are exporters, columns are importers, 3rd dimension is sectors
X=zeros(N,N,K);
for k=1:K
    X(:,:,k) = X_ink(1+(k-1)*N:k*N,:);
end
X_ink=X;

%Other data matrices
V_ik = reshape(V_ik, [N 1 K]);
V_i=sum(V_ik,3);
D_i =reshape(D_i, [N 1 1]);
expenditures_ik = reshape(expenditures_ik,[N,1,K]);
expenditures_nk=permute(expenditures_ik,[2 1 3]);

%Create trade shares, expenditure shares
expenditures_nnk=repmat(expenditures_nk,[N 1 1]);
lambda_ink = X_ink./expenditures_nnk; %trade shares

expenditures_n=sum(expenditures_nk,3);
beta_nk=expenditures_nk./repmat(expenditures_n,[1 1 K]);

%Create sectoral domestic trade shares
dom_lambda_ik=zeros(N,K);
for i=1:N
    for k=1:k
        dom_lambda_ik(i,k) = lambda_ink(i,i,k);
    end
end
        
%Create total sectoral exports
exports_ik=zeros(N,K);
for i=1:N
    for k=1:k
        exports_ik(i,k) = V_ik(i,k)- X_ink(i,i,k);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Gains from Industrial Policy


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Het TE, Het SE (service 0 SE)
theta=h_theta; 
gamma = h_gamma;
counter=0;

gfoip_only_het_0=zeros(N,1);
subsidies_het_0=zeros(N,K);

for i=1:N
    %Defining data vectors and parameters
    V = V_i(i); %total value added
    X=exports_ik(i,:)'/V; %sectoral exports (divided by value added)
    D= D_i(i)/V; %Aggregate trade balance (divided by value added)
    beta = squeeze(beta_nk(1,i,:)); %sectoral expenditure shares
    lambda= dom_lambda_ik(i,:)';  %own expenditure shares by sector
    V_k = V_ik(i,:)'/V; %value added by sector (divided by value added)
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Optimal industrial policy
    %Specifying initial taxes and subsidies
    t_k=0;
    options=optimset('TolX',.005,'TolFun',.0025,'MaxFunEvals',500)
    %Solving the model
    fun=@(s_k)welfare_solver_ces(V_k,D,X,lambda,beta,gamma,theta,rho,t_k,s_k,K);
    %[x fval]=fminsearch(fun,s_k_init,options)
    s_k_init=gamma(2:K)/3;
    [x fval]=fminunc(fun,s_k_init,options)
    gfoip_only_het_0(i)=-fval; 
    subsidies_het(i,:)=[0;x]';
end

%Saving the gains
gains=[countryid gfoip_only_het_0];
csvwrite('gains_20_cen_ces.csv',gains);

